Quantum transport of strongly interacting photons in a one-dimensional nonlinear 

waveguide 



o 
o 

(N 
> 

o 
in 

S 1 
a,: 
-i— > 

cr 



> 



O 



X 



Mohammad Hafezi Q, 1 Darrick E. Chang, 2 Vladimir Gritsev, 3 Eugene Demler, 1 and Mikhail D. Lukin 1 

1 Physics Department, Harvard University, Cambridge, MA - 02138 
2 Center for the Physics of Information and Institute for Quantum Information, 
California Institute of Technology, Pasadena, CA 91125 
3 Physics Department, University of Fribourg, Chemin du Musee 3, 1700 Fribourg, Switzerland 

(Dated: November 25, 2009) 

We present a theoretical technique for solving the quantum transport problem of a few photons 
through a one-dimensional, strongly nonlinear waveguide. We specifically consider the situation 
where the evolution of the optical field is governed by the quantum nonlinear Schrodinger equation 
(NLSE). Although this kind of nonlinearity is quite general, we focus on a realistic implementation 
involving cold atoms loaded in a hollow-core optical fiber, where the atomic system provides a 
tunable nonlinearity that can be large even at a single-photon level. In particular, we show that 
when the interaction between photons is effectively repulsive, the transmission of multi-photon 
components of the field is suppressed. This leads to anti-bunching of the transmitted light and 
indicates that the system acts as a single-photon switch. On the other hand, in the case of attractive 
interaction, the system can exhibit either anti-bunching or bunching, which is in stark contrast to 
semiclassical calculations. We show that the bunching behavior is related to the resonant excitation 
of bound states of photons inside the system. 

PACS numbers: 42.65.-k,05.60.Gg,42.50.-p 



I. INTRODUCTION 

Physical systems that enable single photons to inter- 
act strongly with each other are extremely valuable for 
many emerging applications. Such systems are expected 
to facilitate the construction of single-photon switches 
and transistors [J, [U, networks for quantum infor- 
mation processing, the realization of strongly correlated 
quantum systems using light [1, 0, Q and the investi- 
gation of novel new many-body physics such as out of 
equilibrium behaviors. One potential approach involves 
the use of high-finesse optical microcavities containing a 
small number of resonant atoms that mediate the inter- 
action between photons pH 0|- Their nonlinear proper- 
ties are relatively straightforward to analyze or simulate 
because they involve very few degrees of freedom (i.e., 
a single optical mode) [1, [{J E(|. Recently, an alter- 
native approach has been suggested, involving the use 
of an ensemble of atoms coupled to propagating pho- 
tons in one-dimensional, tightly-confining optical waveg- 
uides [ill, E2> EH • Here, the nonlinearities are enhanced 
due to the transverse confinement of photons near the 
diffraction limit and the subsequent increase in the atom- 
photon interaction strength. The propagation of an op- 
tical field inside such a nonlinear medium (e.g., systems 
obeying the quantum nonlinear Schrodinger equation) is 
expected to yield much richer effects than the case of an 
optical cavity due to the large number of spatial degrees 
of freedom available. Simultaneously, however, these de- 
grees of freedom make analysis much more difficult and 
in part cause these systems to remain relatively unex- 
plored @, EH EE> Efl 03 • We show that the multi-mode, 
quantum nature of the system plays an important role 
and results in phenomena that have no analogue in either 



single-mode cavities or classical nonlinear optics. It is in- 
teresting to note that similar low-dimensional, strongly 
interacting condensed matter systems are an active area 
of research, but most of this work is focused on closed 
systems close to the ground state or in thermal equilib- 
rium [II, El, 53, 5H EH- On the other hand, as will be 
seen here, the relevant regime for photons often involves 
open systems and driven dynamics. 

In this article, we develop a technique to study the 
quantum transport of a few photons inside a finite-length, 
strongly nonlinear waveguide where the light propaga- 
tion is governed by the quantum nonlinear Schrodinger 
equation (NLSE), and apply this technique to study the 
operation of this system as a single-photon switch. In 
particular, we study the transmission and reflection prop- 
erties of multi-photon fields from the system as well as 
higher-order correlation functions of these fields. We find 
that these correlations not only reflect the switching be- 
havior, but reveal some aspects of the rich structure as- 
sociated with the spatial degrees of freedom inside the 
system, which allow photons to "organize" themselves. 
In the regime where an effectively repulsive interaction 
between photons is achieved, anti-bunching in the trans- 
mitted field is observed because of the switching effect, 
and is further reinforced by the tendency of photons to 
repel each other. In the attractive regime, either anti- 
bunching (due to switching) or bunching can occur. We 
show that the latter phenomenon is a clear signature of 
the creation of photonic bound states in the medium. Al- 
though we focus on a particular realization involving the 
propagation of light, our conclusions on quantum trans- 
port properties are quite general and valid for any bosonic 
system obeying the NLSE. 

This article is organized as follows. In Sec. HH we de- 
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scribe an atomic system whose interactions with an opti- 
cal field can be manipulated using quantum optical tech- 
niques such that the light propagation obeys the quan- 
tum NLSE. This method relies upon electromagnetically 
induced transparency (EIT) to achieve resonantly en- 
hanced optical nonlinearities with low propagation losses 
and the trapping of stationary light pulses using spatially 
modulated control fields. Before treating the nonlinear 
properties of the system, we first consider the linear case 
in Sec. IIII1 where it is shown that the light trapping 
technique leads to a field build-up inside the medium 
and a set of discrete transmission resonances, much like 
an optical cavity. In Sec. IIV1 we then investigate the 
nonlinear transport properties of the system such as re- 
flectivity and transmittivity in the semi-classical limit, 
where the NLSE is treated as a simple complex differ- 
ential equation. Here we find that the presence of the 
nonlinearity causes the transmission resonances to shift 
in an intensity-dependent way - the system behaves as 
a low-power, nonlinear optical switch, whose behavior 
does not depend on the sign of the nonlinear interac- 
tion. In Sec. El we present a full quantum formalism 
to treat the NLSE transport problem in the few-photon 
limit. Sec. I VII is dedicated to analytical solutions of the 
NLSE with open boundary conditions when the system is 
not driven. In particular, wc generalize the Bethe ansatz 
technique to find the resonant modes of the system, which 
help to elucidate the dynamics in the case of the driven 
system. The driven system is studied in Sec. IVII1 where 
numerical solutions are presented along with a detailed 
study of the different regimes of behavior. In particular, 
we find that the correlation functions for the transmitted 
light do depend on the sign of the nonlinear interaction, 
in contrast to what the semi-classical calculations would 
suggest. We conclude in Sec. IVIII1 



II. MODEL: PHOTONIC NLSE IN ID 
WAVEGUIDE 

In this section, we consider the propagation of light 
inside an finite-length atomic medium under EIT con- 
ditions and with a Kerr nonlinearity. We also describe 
a technique that allows for these pulses of light to be 
trapped within the medium using an effective Bragg 
grating formed by additional counter-propagating optical 
control fields. We show that in the limit of large optical 
depth the evolution of the system can be described by a 
nonlinear Schrodinger equation. 

Following Ref. [6(, we consider an ensemble of atoms 
with the four-level internal structure shown in Fig. [1] 
which interact with counter-propagating quantum fields 
with slowly- varying envelopes £± inside an optical waveg- 
uide. These fields are coupled to a spin coherence 
between states \a) and |c) via two classical, counter- 
propagating control fields with Rabi frequencies f2± 
largely detuned from the \b) — > |c) transition. The case 
where the fields propagate only in one direction (say in 



the "+" direction) and where the detuning is zero cor- 
responds to the usual EIT system, where the atomic 
medium becomes transparent to £ + and the group veloc- 
ity can be dramatically slowed due to coupling between 
the light and spin wave (so-called "dark-state polari- 
tons" ) [23j . On the other hand, the presence of counter- 
propagating control fields creates an effective Bragg grat- 
ing that causes the fields £± to scatter into each other. 
This can modify the photonic density of states and cre- 
ate a bandgap for the quantum fields. This photonic 
bandgap prevents a pulse of light from propagating and 
can be used to effectively trap the light inside the waveg- 
uide [U I25I ]. The trapping phenomenon is crucial be- 
cause it increases the time over which photons can inter- 
act inside the medium. The presence of an additional, 
far-detuned transition |c) — > \d) that is coupled to £± 
leads to an intensity-dependent energy shift of level |c), 
which translates into a Kerr-type optical nonlinearity 

We now derive the evolution equations for the quan- 
tum fields. We assume that all atoms are initially in their 
ground states \a). To describe the quantum properties 
of the atomic polarization, we define collective, slowly- 
varying atomic operators, averaged over small but macro- 
scopic volumes containing N z ^> I particles at position 
z, 

1 Nz 

z i— 1 

The collective atomic operators obey the following com- 
mutation relations, 



[a a p(z),a flI/ (z')} = —8(z - z')(8 0ll a av (z) - S^a^iz)), 

(2) 

while the forward and backward quantized probe fields in 
the z direction obey bosonic commutation relations (at 
equal time), 

[£\(z),£{(z')} = 5(z-z'). (3) 

The Hamiltonian for this system in the rotating frame 
can be written as 



H = - - J [a bc (z)(n + e^ z + n^e- tk ^) + H.c.} (4) 

+ gV2^[a ba {z)(£+e lkoZ + ^ e - lk ° z ) + H.c] 
+ gV2^[a dc (z)(£+e tk '-> z + £_ e - lk ° z ) + H.c] 
+ Ai& bb (z) + A 3 a cc (z) + (A 2 + A 3 )add(z)dz 

where g = l^^j^^j is the atom-field coupling strength, 

/i is the atomic dipole matrix element, and A is the effec- 
tive area of the waveguide modes. For simplicity, we have 
assumed that the transitions a-b and c-d have identical 



coupling strengths g and have ignored transverse varia- 
tion in the fields. The terms Aj denote the light field- 
atomic transition detunings as shown in Fig. [lja). fc c is 
the wavevector of the control fields, while fco — nbto a b/c 
characterizes the fast-varying component of the quan- 
tum field and rib is the background refractive index. We 
also define no = N/L as the linear density of atoms in 
the z direction, and v g = cfl 2 /2irg 2 no as the group ve- 
locity that the quantum fields would have if they were 
not trapped by the Bragg grating (we will specifically 
be interested in the situation where Q + = f2_ = f2). 
Following Ref. [23[, we can define dark-state polariton 
operators that describe the collective excitation of field 
and atomic spin wave, which in the limit of slow group 
velocity r) = ^f- > 1 are given by #± = gv ^"° ^±- 
These operators obey bosonic commutation relations 
[^±(z)A±(z')} = S(z - z'). The definition of the po- 
lariton operators specifies that the photon flux entering 
the system at its boundary is equal to the rate that po- 
laritons are created at the boundary inside the system - 
i.e., c(£+£+) = v g In other words, excitations 

enter (and leave) the system as photons with velocity 
c, but inside the waveguide they are immediately con- 
verted into polariton excitations with group velocity v g . 
Field correlations will also be mapped in a similar fash- 
ion - in particular, correlation functions that we calcu- 
late for polaritons at the end of the waveguide z — L will 
also hold for the photons transmitted from the system. 
The total number of polaritons in the system is given 
by Afpoi. = f + The optical fields 

coupled to the atomic coherences of both the a — b and 
c—d transitions are governed by Maxwell-Bloch evolution 
equations, 



d d \ 

— ± c— J £t (z, t) = igV2Trn (a- ab + a cd 



)• (5) 



Similar to the photonic operators, the atomic coher- 
ences can also be written in terms of slowly- varying com- 
ponents, 



<Jab — <J „h e T cr , e , 



a+e ikoZ + a-,e~ ik ° z . 



(6) 
(7) 



We note that higher spatial orders of the coherence are 
thus neglected. In practice, these higher orders are de- 
stroyed due to atomic motion and collisions as atoms 
travel distances greater than an optical wavelength dur- 
ing the typical time of the experiment 27]. Alternatively, 
one can use dual-V atomic systems that do not require 
this approximation [281 ]. 

In the weak excitation limit (a aa — 1), the population 
in the excited state \b) can be neglected, (abb)~0- In this 
limit, the evolution of the atomic coherence is given by 



+ igV2^£± + ifl±a ac e ±iAkz 



(8) 
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FIG. 1: (a) Four-level atomic system for creating strong non- 
linearity. Counter-propagating control fields modulate the 
EIT for the forward- and backward-propagating probe, and 
c) — » \d) transition gives rise to a Kerr- type nonlinearity. (b) 
The light is confined in the transverse direction due to the 
presence of the waveguide and experiences an effective Bragg 
grating due to the presence of the counter-propagating light. 



where Afc — k c — fco and T is the total spontaneous emis- 
sion rate of state b (for simplicity we also assume that 
state d has an equal spontaneous emission rate). In the 
adiabatic limit where g^/2n(a ac ^±) *C T, the coherence 
a a d can be approximated by 



a ad ~ / V2 7 ac , r (* + e+*°' + *_ e -*>*). (9) 
-A 2 - A 3 - 

Therefore, the spin wave evolution can be written as, 

& ac = iA 3 & ac + i{cr+n* + e- iAkz + *- b n*_e lAkz ) 
2irig 2 



L {£{& + + £{&-)a ac (10) 

1 2 

We now consider the situation where f2 + = f)_ = £7, 
such that the counter-propagating control fields form a 
standing wave. In the adiabatic limit [231 ] , and keeping all 
terms up to third order in the quantum fields, substitut- 
ing these results into Eq. ([5|) and simplifying yields the 
following evolution equations for the dark-state polariton 
operators, 

(cd z + d t )*+ = -|(* + -#_)_|a t (# + + #_) 
+ (* t + + * t _)(*+ + *_)*+] (ii) 
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(-cd 2 +9 t )*- 



+ (* f + + * t _)(*+ + #_)#_], (12) 
where the linear dispersion is characterized by £ = 

27r 2 n 

-iAi+r/2 ' r ^^ ie nonunear ity coefhcient is given by the 

2 

single photon AC-Stark shift: A„ = jfA^+TrJTj- We no * e 
that the wave-vector mismatch Afc has been compen- 
sated for by a small extra two-photon detuning equal 
to (-Akc/rj). 

The above equations describe the evolution of two cou- 
pled modes. It is convenient to re-write these equations 
in terms of the anti-symmetric and symmetric combina- 
tions A = (#+ - #_)/V2 and 5 = (*+ + *_)/v / 2. For 
large optical depths, we then find that the anti-symmetric 
mode adiabatically follows the symmetric mode, A ~ 
— (c/£)d z S. In this limit, the evolution of the whole sys- 
tem can be described by a single nonlinear Schrodinger 
equation, 



d q 



£ dz 2 



S + 8iA n S^ S 2 = 0. 



Physically, the coupling between ^± induced by the 
Bragg grating causes them to no longer behave in- 
dependently, much like the two counter-propagating 
components of an optical cavity mode. We can 
write the above equation in dimcnsionlcss units by 
introducing a characteristic length scale L co h 
c/|Im[£]| = c(A? + r 2 /4)/27rg 2 n |A 1 | and time scale 
t coh = 77/|Im[e]| = (A 2 + r 2 /4)/2n 2 |A!|. L coh corre- 
sponds to the length over which the field acquires a tt- 
phase in the propagation. The dimensionless NLSE then 
reads 



.as 



i d 2 s 

2m dz 2 



2/tS t S 2 , 



(14) 



where for Ai < 0, the effective mass is to = |(1 + » 2 |Ai| ) 
and the nonlinearity coefficient is k — 2 -^ L — 7rg ^ c 



A 2 +ir/2- 



Note that ^±(z, t) and S(z, t) are also in units of y L co ) v 

such that [S+(z), S — 5(z — z'). For simplicity, 
we omit tilde superscripts in the following. We can 
also write the nonlinear coefficient as k 



4(A 2 +ir/2) ' 

where we have identified Tm = 4:irg 2 /c as the spon- 
taneous emission rate into the guided modes (ri£><r). 
We are primarily interested in the limit | Ai : 2 |3>T such 
that to, k are mostly real and the evolution is dispersive. 
Note that in this notation, the anti-symmetric combi- 
nation of forward and backward polaritons is given by 
A~ -i/2md z S~ -id z S. 



III. LINEAR CASE: STATIONARY LIGHT 
ENHANCEMENT 

In this section, we investigate the linear transmission 
properties of the signal field as a function of its frequency. 
The control field leads to a Bragg grating that couples 
the forward and backward components of the signal field 
together. We show that the system therefore acts as an 
effective cavity whose finesse is determined by the optical 
density of the atomic medium. 

For the linear case (k — 0), it is sufficient to treat 
the forward and backward field operators as two complex 
numbers. In the slow light regime (77 3> 1), the coupled 
mode equations (Eqs. I12[) can be written in the Fourier 
domain, with our dimcnsionlcss units, as 



-5($+ + $_) + im($ + -$_) (15) 



-«9 Z $_ = + $_) - im($+ - $_), (16) 

where *+(z,t) = $ + (z, 5)e~ i&T and t) = 

<&_(;z, <5)e~ ll5T and 5 is the dimensionless two-photon de- 
tuning 5 — A 3 t coh . We specify that a classical field 
<f> + (z = 0,5) = a enters the system at z — with 
(13) no input at the other end of the system (z = d), 
$_(z = d,S) = 0, as shown in Fig. QJb). We note 
that d — L I is the length of the system in units of 
the coherence length introduced earlier. For negligible 
losses (|Ai|3>r ) and Ai < 0, to ~ 1/2 and the profile of 
forward-going polaritons inside the system will look like: 

L $+(z,8) 2iVS cos[{d - z)VS] + (1 + S) sin[(d - z)VS] 



a 2i\f5 cos{d\/5] + (1 + 5) sm[dy/S\ 

(17) 

Therefore, for a system with fixed length d, the trans- 
mission coefficient varies with the frequency of the in- 
cident field, with transmission resonances occurring at 
the values \fSod = n7T ( n is an integer). At these 
resonances, the system transmittance is equal to one 
(\3>(d,5)\ = |$(0,<5)|) and a field build-up occurs in- 
side the medium with a bell-shaped profile, similar to 
a cavity mode (see Fig. The positions of these reso- 
nances (quadratic in n) reflect the quadratic dispersion 
in Eq. (TT4")) . Note that in real units, the positions of the 
resonances will depend on the amplitude of the control 
field, since A3 = <5^j^p-- In the limit of a coherent opti- 
cally large system (d> 1), the intensity amplification in 
the middle of the system is equal to (d/27r) 2 for the first 
resonance. In other words, the Bragg scattering creates a 
cavity with an effective finesse proportional to the square 
of the coherent length of the system [T oc d 2 ). 

We now derive the width of the first transmission reso- 
nance. For small variations <5o ± 5b around the resonance 
frequency, we can write 



$+(<*) 
*+(o) 



= -1- 



45, 



3/2 



S b + 77^ 2 + 0(6%). (18) 
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FIG. 2: Linear case: (a) Transmittivtiy as a function of two- 
photon detuning. Transmission peaks are attenuated because 
of linear loss on a) — > \b) transition which are plotted for 
three different loss rates f3 = d T/Ai. (b) When the system 
is tuned on a transmission resonance (\f~5od — nn), the field 
inside the medium is amplified. 



Therefore, the width of the resonances (say where it 
drops by half) is given by 



25 b 



5 3/2 



<5> 



(19) 



We have kept terms up to second order in 5b, since the 
first order term does not give a decreasing correction to 
the transmittancc. While we have previously ignored ab- 
sorption (as determined by the real part of £), we can 
estimate that its effect is to attenuate the probe beam 
transmission by a factor f3 = d T/|Ai|. As it is shown in 
Fig. [3J for large optical densities, (3 can fully character- 
ize the transmission coefficient on resonance (5 = 5 ). In 
particular, for a fixed [3, the resonant transmission is con- 
stant for any large optical density. In other words, since 
the optical depth of the system is given by OD = rf^M, 
the transmittivity of the system remains constant for any 
OD with the choice |Ai| = Ty/OD/p. In this case the ef- 
fective cavity finesse for the system becomes proportional 
to the optical density, i.e., T cx OD. 

The total number of polaritons in the system can be 
estimated by, 



Kol = [ d \<P+(z)\ 2 + \<P-(z)\ 2 dz 

Jo 



(20) 



(d 2 +7T 2 ) 2 , , d 3 , , 

4^l* + (0)|^^|$ + (0)| 2 . 



This again shows that the polaritons experience many 
round trips inside the system before exiting. In par- 
ticular, if we define the average intensity inside the 
medium as |<i>^ ue | 2 = Af po i/d, then we readily observe 
that the intensity of the polariton field is amplified 
inside the medium by the square of the system size 
(|$^ e | 2 /|$ + (0)| 2 = (d/2n) 2 ) - i.e. the finesse is pro- 
portional to the optical density (OD). 

The original proposal for observing an enhanced Kerr 
nonlinearity with a four-level atomic system using EIT 




200 400 600 800 1000 
OD 

FIG. 3: For large optical densities, the transmission on reso- 
nance (5 = So), only depends on j3 = OD • 



makes use of an optical cavity to enhance the nonlinear- 
ity [![. However, as pointed out in Rcf. [9], the scheme 
suffers from some inaccuracies in the effective Hamilto- 
nian. More specifically, in Ref. Q, the effective Hamil- 
tonian was evaluated at the center of the EIT trans- 
parency window. However, in practice, EIT dramatically 
decreases the cavity linewidth because of the large dis- 
persion that accompanies the vanishing absorption [29| ; 
this causes photons at frequencies slightly shifted from 
the central frequency to be switched out of the cavity. 
This leads to an extremely small allowable bandwidth 
for the incoming photons [9J] and was neglected in the 
original analysis. We emphasize that the analysis pre- 
sented here takes into account the dispersive properties 
of the medium, as we have included the field dynamics 
up to second order in the detuning from resonance (this 
accounts for the effective mass of the photons in our sys- 
tem). We verify the consistency of this derivation in Ap- 
pendix [X] by solving the linear system including full sus- 
ceptibilities. It is shown that the results are consistent 
near the two-photon resonance (i.e., frequencies around 
5 = 0). 



IV. SEMI-CLASSICAL NONLINEAR CASE 
A. Dispersive regime 

In this section, in contrast to the previous section, we 
include the nonlinear term in the evolution equations to 
investigate its effect in the semi-classical limit (where the 
fields are still treated as complex numbers). In this pic- 
ture, the effect of nonlinearity causes the transmission 
peaks to shift in frequency in an intensity-dependent way 
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to the left or right depending on the sign of the nonlin- 
earity coefficient k. We show that when |k|0£> ^> 1, the 
magnitude of the shift is large even at intensities corre- 
sponding to that of a single photon. In this regime, we 
expect that the system can act as a single-photon switch 
and that signatures of quantum transport will become 
apparent (the quantum treatment is described in Sec fV)) . 

Because of the self-phase modulation term in the evo- 
lution equations (Eqs. I12p , the forward and backward 
fields acquire a phase shift proportional to their inten- 
sity. Moreover, due to the conjugate-phase modulation 
terms, each field undergo an extra phase shift propor- 
tional the intensity of the other field. Classically, this 
yields a frequency shift in the transmission spectrum 
when the nonlinearity is small. The shift in the trans- 
mission peak can be approximated by AS ~ 2ft|<I>™ e | 2 

where |$™ e | 2 ~ ^|$ + (0)| 2 is the average intensity of 
polaritons in the system. Suppose that we want the non- 
linearity to be strong enough to shift the transmission 
peaks at least by half of their widths, AS ~ |^ 3 ^ 2 - Then, 
from Eq. tflU)) this condition can be written as 



i $c +w = n) 5 o 



(21) 



On the other hand, according to Eq. (f21j) . we can write 
this condition in terms of the critical number of polari- 
tons inside the system, 



7V cr , = —— 



(22) 



Since the nonlinearity coefficient is given by the light 
shift on the |c) — > \d) transition, in the dispersive 
regime (A 2 T), we have k — Thus, we ex- 

pect to have substantial nonlinearities at the level of one 
polariton (i.e., one incoming photon), M"\ = 1, if 



r a. 



ID 



(23) 



where Tm is the rate of spontaneous emission rate into 
the guided modes. Strictly speaking, note that a single 
photon cannot actually have a nonlinear phase shift (as 
correctly derived later using a fully quantum picture); 
however, we can still use the results of this semiclassi- 
cal calculation to qualitatively understand the relevant 
physics. 

We can also rewrite the above condition in term of the 
optical density [OD = d4r) needed in the system. From 
the linear case, we know that an optimal detuning, for 
a transmission of 90%, should satisfy d-£- ~ 0.5. Then, 
Eq. (f23|) can be written as 



OD = 27r d 



r A 2 
r 1D r 



(24) 



Taking for example a system where A 2 ~ 5r and ^4P- ~ 
0.1, nonlinearities at a few-photon level can be observed 
for an optical density OD ~ 6200. 

First, let us consider the case of positive k. In Fig. 
SI we observe that at large enough optical density, the 
system can have very different transmission spectra for 
low and high intensities that classically correspond to 
having one and two polaritons (photons) inside the sys- 
tem, respectively. Although we have ignored the quan- 
tization of photons in this section, we can develop some 
insight into the transmission properties of one- and two- 
photon states. Loosely speaking, if we fix the input field 
frequency to lie at the one-photon (linear) transmission 
peak (So), the system would block the transmission of 
incident two-photon states. More realistically, suppose 
we drive the system with a weak classical field (coher- 
ent state), which can be well-approximated as containing 
only zero, one, and two-photon components. We then 
expect that the one-photon component will be transmit- 
ted through the system, while the two-photon component 
will be reflected, leading to anti-bunching of the trans- 
mitted light. We note that the general spirit of this con- 
clusion is sound; however, the correct description of the 
system is achieved by taking into account the quantiza- 
tion of photons which is presented in the next sections. 

A similar analysis holds for the case of negative k. Note 
that the sign of k depends on the detuning of the photonic 
field from the atomic transition |c) — > \d), which can eas- 
ily be adjusted in an experiment. This is in contrast to 
conventional nonlinear optical fibers and nonlinear crys- 
tals, where the nonlinearity coefficient is fixed both in 
magnitude and sign. We find that a negative nonlinearity 
simply shifts the transmission spectrum in the opposite 
direction as for the positive shown in Fig. [5l but 

all other conclusions remain the same. In particular, we 
would expect anti-bunching to occur for this case as well, 
when a weak coherent field is incident with its frequency 
fixed to the linear transmission resonance. Surprisingly, 
the quantum treatment fSec lVII[) . shows that the above 
conclusion is wrong and system behaves very differently 
for negative nonlinearity. We show that this difference in 
behavior can be attributed to the presence of additional 
eigenstates (photonic bound states) in the medium and 
their excitation by the incident field. 

For even larger nonlinearities or intensities, the trans- 
mission spectrum can become even more skewed and ex- 
hibit bistable behavior, as similarly found in Ref. [13] 
in the context of transport of Bose-Einstein condensates 
in one dimension. There, the classical NLSE (Gross- 
Pitaevskii equation) was solved to find the mean-field 
transport properties of a condensate scattering off a po- 
tential barrier. 

Instead of considering the switching effect as a func- 
tion of number of photons inside the medium, we can 
also consider the number of photons that need to be sent 
into the system. Clearly, to have a well-defined trans- 
mission amplitude without substantial pulse distortion, 
the incident pulse must be long enough so that it fits 
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within the bandwidth of the system resonance, as given 
in Eq. (fT!)]) . To be specific, we consider an input pulse 
whose duration is equal to the inverse of the bandwidth, 
tb = (§) 3 t co h- We can relate the number of incoming 
photons to an average incident intensity: 



|<M0)| 2 =AU 



tb 



coh 



JVpol. 



(25) 



Now, since the number of incident photons and incom- 
ing polaritons are the same, we can assign an average am- 
plitude to any incoming photons number by Eq. (|25D , and 
evaluate the transmission. Hence, we can evaluate the 
number of incident photons needed to observe a signifi- 
cant nonlinearity and saturate the system. Fig. [5] shows 
the transmittivity of the nonlinear system as a function 
of number of photons in the incoming wavepacket. We 
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observe that for high optical densities (OD > 1000), the 
transmittivity drops as the number of incoming photons 
increases and the system gets saturated for even few pho- 
tons. 



B. Dissipative Regime 

In this section, we investigate the system in the pres- 
ence of nonlinear absorption, where k is imaginary. The 
nonlinear dispersion of the previous case can simply be 
turned into nonlinear absorption by setting the nonlinear 
detuning to zero (A2 = 0, k = -gif^)- I n the quantum 
picture, this term does not affect the one-photon state, 
while two-photon states can be absorbed by experienc- 
ing three atomic transitions, \a) — > \b) — > |c — > \d), and 
subsequently being scattered from excited state \d) [3ll |. 
We consider the quantum treatment of absorption later 
and first study the semiclassical limit here. 

The presence of nonlinear absorption suppresses the 
transmission of multi-photon states through the medium 
by causing them to decay. This suppression becomes 
stronger for higher intensities as shown in Fig. [71 We 
have used the same optical density (OD) and ID confine- 
ment (Tin /r) as in Fig. |4j We observe that the effects of 
nonlinear absorption are stronger than that of nonlinear 
dispersion studied in Sec lIV Al since it occurs at reso- 
nance (A2 = 0) where the atomic response is strongest. 
It is thus possible to observe its effect at even lower inten- 
sities, corresponding to effective photon numbers two or- 
ders of magnitude smaller than the dispersive case. Much 
like the dispersive case, the suppression of transmission 
of multi-photon components should yield anti-bunching 
in the transmitted field. In this case, however, these com- 
ponents are simply lost from the system (as opposed to 
showing up as a bunched reflected field). 
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the field number operator S'S, manifolds with different 
field quanta are decoupled from each other inside the 
medium. Therefore, the evolution for the one-photon 
and two-photon manifolds under the NLSE Hamiltonian 
can be written as, 



■ 9 M 



i ( d 2 a 2 



2m \dzf dz 2 
+ 2n<j)(z 1 , z 2 ,t)S(zi - 



Zl,Z 2 ,t) 

*a) (27) 
(28) 



However, the system is driven with an input field at 
z = 0, which allows different manifolds to be coupled 
at the boundaries. This is analogous to fiber soliton 
experiments where a classical input field mixes quantum 
solitons with different photon numbers [13, [H, [34]] • In 
particular, for a classical input field, 



V. QUANTUM NONLINEAR FORMALISM: 
FEW-PHOTON LIMIT 

In this section, we describe a quantum mechanical ap- 
proach that enables one to solve the problem of quantum 
transport of a small number of photons through the finite, 
nonlinear system described in Sec|lTJ This few-photon 
number limit is of particular interest since it captures 
the physics of single-photon switching. 

We find it convenient to study the dynamics of the sys- 
tem of photons in the Schrodinger picture, where one can 
explicitly solve for the few-body wave functions. This ap- 
proach is made possible by truncating the Hilbert space 
so that only subspaces with n max photons are less are 
present. In the following, we will consider the case where 
n ma x = 2, although our analysis can be easily extended 
to cover any other value. This truncation is justified 
when the incident coherent field is sufficiently weak that 
the average photon number is much smaller than one in- 
side the system (|q;o | 2 cZ 3 <C 1, where «o is the amplitude 
of the incoming field). Thus, we can write the general 
state of the system as: 



* I dz 1 dz 2( f>(z 1 ,z 2 ,t)S\z 1 )S\z 2 )\0) (26) 
dzd(z,t)S^(z)\0) + e\0). 



The first, the second and the third term correspond to 
two-photon, one-photon and vacuum state, respectively. 
Note that because of bosonic symmctrization, 4>(z\, z 2l t) 
should be symmetric in z\ and z 2 . This formalism allows 
us to capture any non-trivial spatial order between 
photons in our system {e.g., the de- localization of two 
photons as represented by the off-diagonal terms in 
c6(zi, z 2 )). Since the NLSE Hamiltonian commutes with 



*+(* = 0)|V(t)) = a(t)\4>(t)) (29) 

*_(2 = d)Mt)) = 0\m), (30) 

which corresponds to a coherent state with (possibly 
time-dependent) amplitude a(t) as an input at z = 0, 
and no input (i.e., vacuum) at z = d. Since we specify 
that the input coherent field is weak (a <C 1), the 
amplitude of the vacuum state is almost equal to one 
(e ~ 1). The annihilation operator in these equations 
reduces the photon number on the left-hand side by 
one. Thus, such boundary conditions relate different 
photon subspaces whose photon number differ by one, 
e.g. the two-photon and one-photon wavefunctions. In 
the adiabatic limit where the anti-symmetric part of the 
field (A = (^ + — ^>-)/\/2) follows the symmetric part 
(S = (* + + *_)/ v / 2), we have 



Therefore the boundary conditions at z — can be re- 
written as 

^ I d Zl dz 2 [5-^ ( 9^] z=0 5 t (z 1 )5 t (z 2 )c/»(zi,Z2,t)|0) 



l 



dz0(z,t)S*\O), 



dz [S-—d z S] z=0 S\z)6(z,t)\0)=a\0). 
2m 



(32) 



Using the identity J[d z S(z), &(z')]f(z')dz' = d z f(z), 
the boundary conditions on the one-photon and two- 
photon wave functions can be written as: 



(zi ^0,z 2 ,t)~^-d^^i,z 2 ,t)\ zl=0 = -j=0(z2,t) 

0(z = 0,t)- ^-d z 9(z = 0,t) = V2a, (33) 
2 m 
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where c^ 1 ' acts on the first parameter. This type of open 
boundary condition is known as a Robin or mixed bound- 
ary condition, which involves a combination of both the 
function and its derivative. In the present case, the open 
boundary conditions allow particles to freely enter and 
leave the system. We emphasize that this process is 
noise-less, in that the loss of population from the interior 
of our system is related by our boundary condition equa- 
tions to the flow of particle current through the system 
boundaries. This is in contrast to an optical cavity, for in- 
stance, where photons inside the cavity leak dissipatively 
into the environment (35| . Similarly the boundary con- 
dition at z = d reads 



<j>(d,z,t) + —dW<j>(d,z,t) 
2m 



9{d,t) 



2m 



d z 6{d,t) 



(34) 
(35) 



Given the boundary conditions and the equations of mo- 
tion in the interior, we can completely solve for the pho- 
ton wavefunctions. 

Once the wavefunctions are determined, it is possible 
to determine the intensity profile as well as any other 
correlation function for the photons. For example, the 
intensity of the forward-going polariton is 

I(z,t) = (V>(t)|$Uz)*+(*)Mt)) (36) 
= (l\¥ + (z)i> + (z)\l) + (2\¥ + (z)^ + (z)\2), 

where \ j) denotes the component of the total wavefunc- 
tion |?/>(i)) containing j photons. The first and second 
terms on the right thus correspond to the one- and two- 
photon contributions to the intensity. By re-writing ex- 
pressions in terms of 5* instead of vfr-), we obtain: 



(z)-—d z 9(z) 
2m 



(37) 



dz' 



{ z ', z )-^-d^<f>{z',z) 
2m 



and physically characterizes the relative probability of 
detecting two consecutive photons at the same position 
z. If this quantity is less (greater) than one, the pho- 
tonic field is anti- bunched (bunched). In particular, if 
g<i{z) — 0, the field is perfectly anti-bunched and there 
is no probability for two photons to overlap in position. 
In our truncated Hilbcrt space, gi(z) of the transmitted 
field is given by 



g 2 (z = d) 



<2|*£(d)*3.(d)|2) 



(l|*Ud)*+(d)|l) + (2|*t(d)*+(d)|2) 



We note that this expression can be simplified, since at 
z = d, we have = -j-(S + ^d z S) = and = 

V2S. Therefore, 



92(d) 



4\(/>(d,d)\' 



(\9(d)\ 2 + Ajdz'\^z',d)\ 2 y 



(41) 



We can also evaluate the stationary two-time correlation, 
which is defined in the Heisenberg picture as: 



92{z, 



(tP\*1(z,0)*1(z, 



where the denominator is simplified in the stationary 
steady-state regime. This correlation function charac- 
terizes the probability of detecting two photons at posi- 
tion z but separated by time r. We can re-write g 2 {z, r) 
in terms of wavefunctions in the Schrodinger picture in 
the following way. We first note that the expression 
|V>(0)) = ^+(z, 0)\ip) appearing in the equation above 
can be thought of as a new wavefunction, which describes 
the state of the system after a photon is initially detected 
at time t = and position z. This new state naturally 
has one less photon than the original state, and by sim- 
plifying the expressions, it can be written as: 



Similarly, the second-order correlation function for the 
forward field is 



[z,z)--d^<j>{z,z) 
m 



(38) 



Am 2 



which in our truncated space only depends on the two- 
photon wave function. Now, we evaluate the normalized 
second-order correlation function g 2 {z), which character- 
izes the photon statistics of an arbitrary field. This func- 
tion takes the form 



92{z) 



\¥ + 2 (z)^l(zM 



(1>\& + (z)$+(z)\1;) 



(39) 



|^(0)) 



r (z')S^(z')\0} + e new \0) (43) 



where the new one-photon and vacuum amplitudes are 
given by 



— =■ = cj>(d,z',t = Q)-— d^cj>(d,z',t = 0) (44) 
V2 2m 



V2 



(6(d,t = 0)-id6(d,t = 0)). (45) 



Here we have assumed that z = d, since we are interested 
in the transmitted field. Now, Eq. (|4"2"|) can be written as 

g 2 (d,r) = (^(0)|*! h (d 1 r)* + ( ( i ) r)|^(0))/K^(0)|^(0))| a . 

(46) 
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The numerator describes the expectation value for the in- 
tensity operator I{t) = ^ + (d, r)^!+(d, r) in the Heisen- 
berg picture given an initial state |^(0)). However, we 
can easily convert this to the Schrodinger picture by mov- 
ing the evolution from the operator to the state, i.e., 
by evolving | -0 (0)) under the same evolution equations 
(Eqs. [27ll28f and boundary conditions (Eq. [33)) that we 
used earlier. Therefore, the correlation function g 2 (z,r) 
will be given by: 

$(t)\¥ + (z)* + (zM(t)) 
92{z,t) = ± (47) 

<^(0)l*t (*)*+(*) 1^(0)) 



VI. ANALYTICAL SOLUTION FOR NLSE 
WITH OPEN BOUNDARIES 



and similarly for z = d, 

0(d) + 7^-d z 8(d) = 0. (50) 
Zm 

We look for stationary solutions of the form 9(z, t) = 
e~' lSt 9(z), where 9(z) = Asin(fcz) + Bcos(kz). For sim- 
plicity, we assume m = 1/2. Therefore, we recover the 
quadratic dispersion relation S = k 2 . The values of k 
are allowed to be complex to reflect the open nature of 
our boundary conditions, which allows particles to freely 
enter or leave. By enforcing the boundary conditions we 
get a set of equations for the coefficients A, B, 

B - iAk = 0, 

(A - iBk) sin(fcd) + (B + iAk) cos(kd) = 0, 

which yields the characteristic equation for finding eigen- 
modes and eigen-energies of system, 



In this section, we show that a NLSE system with 
open boundary conditions yields analytical solutions in 
absence of an outside driving source (ag = 0) . To obtain 
the analytical solutions, we use the Bethe ansatz tech- 
nique [H, HI] . This ansatz specifies that the eigenstates 
consist of a superposition of states in which colliding par- 
ticles exchange their wavenumbers ki. Unlike the typical 
formulation, the values of fcj here can be complex to re- 
flect the open nature of our boundary conditions, which 
allow particles to freely enter or leave. In particular, 
we present the one-, two- and many-body eigenmodes of 
the system along with their energy spectra. Finding cer- 
tain eigenmodes of the system (e.g., bound states) helps 
us understand the correlation functions and also spatial 
wavefunctions which are numerically calculated later in 
Sec lVIII for a driven system. 



A. One-particle problem 

First, we calculate the fundamental modes for the one- 
particle states. These modes are of particular interest 
when we later want to construct the many-body wave- 
function of the interacting system in the absence of an 
input field. 

Specifically, we want to find solutions of the 
Schrodinger equation for a single particle in a system of 
length d, 



■ 9 nl , 1 d 2 „. , 



dt 



(48) 



subject to open boundary conditions. The boundary con- 
dition for the undriven system at z = is given by 



9(0) 



2m 



d z 9(0) = 



(49) 



2ikd 



1 



(51) 



Therefore the normalized corresponding wave function 
for each allowed k will be: 



(z) = A(ain(kz) + ikcos(kz)), 



(52) 



A = 



4k 



2dk(l + k 2 ) + (k 2 - 1) sin(2dfc) 



We note that in the limit of large optical density d 3> 1, 
the lowest energy modes of the open system are very 
close to those of a system with closed boundary condi- 
tions, whose characteristic equation is given by kd = nn. 
For example, at d = 100, the wave number correspond- 
ing to lowest energy is fc~0.0314-i0.00063~ tt/100. We 
note that the many-body solutions of the system in the 
presence of very strong interactions (large k) can be con- 
structed from these single-particle solutions and proper 
symmetrization, as we show in Sec lVI El 



B. Two- particle problem 

In this section, we study the problem of two particles 
obeying the NLSE with mixed boundary conditions. We 
wish to solve 



1_ fd_ 2 

2m \dz 2 
2k0(zi, z 2 )S(zi 



dz 2 



Z2), 



(53) 



where E is the energy of the system and can be complex. 
Again, we assume the mass is entirely real, m = 1/2. 

We should note that the conventional method of sep- 
aration of variables cannot be applied in this case. The 
reason for this can be understood in the following way. 
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On one hand, if we ignore the delta interaction term 
in the evolution equation of the two particles, finding 
the eigenfunctions is essentially equivalent to solving the 
Laplace equation in a box with mixed boundary condi- 
tions. Therefore, for this problem the natural separation 
of variables involves solutions given by products of func- 
tions f(zi) and g(z 2 ). On the other hand, if we neglect 
the boundaries, the problem of two particles interacting 
at short range can be solved by utilizing the center of 
mass and relative coordinates and invoking solutions in- 
volving products of functions f(z± + Z2) and g(z\ — z-i). 
We immediately see that the two sets of solutions are ir- 
reconcilable and thus separation of variables is not appli- 
cable when both the boundary conditions and interaction 
term are present. 

We thus take a different approach, using a method 
similar to the Bethe ansatz method for continuous, one- 
dimensional systems [l8j]. Specifically, we solve the 
Schrodinger equation in the triangular region where < 
z\ < Z2 < d, and we treat the interaction as a boundary 
condition at Z\ = z-i- In other words, when two particles 
collide with each other at Z\ = Z2, they can exchange mo- 
menta, which is manifested as a cusp in the wave function 
at z\ = Z2- Hence, for the boundary conditions in this 
triangular region, we have 



0(O,z 2 )-i9 Zl 0(O,z 2 ) = (54) 
<p(zi, d) + id Z2 cp(zi, d) = (55) 
(d Z2 - d Zl )<j>(z 1 ,z 2 )\ Z2=Zl = K(f>(zi,z 2 )\ Z2=Zl (56) 

We note that the last boundary condition is deduced from 
integrating Eq. (|53[) across z\ = Z2 and enforcing that the 
wavefunction is symmetric, 

(d Z2 -d Zl )<p(zi,z 2 )\ Z2=z + - {d Z2 -d Zl )(j>(zi,z 2 )\ Z2 = z - 

= 2K<f)(zi,z 2 )\ Zl=Z2 . (57) 

Inside the triangle, the solution consists of superpositions 
of free particles with complex momenta. Since particles 
can exchange momenta when they collide at Z\ = Z2, we 
should consider solutions of the following form, 



<j>(zi,z 2 ) =^A e e ieiklZl+le2k2Z2 +B e e ieik2Zl+te2klZ2 

(58) 

where the summation should be performed on all sets of 
signs e = ±1. Given the terms containing A e , the terms 
B € then arise from the scattering of the particles off each 
other. Let's first consider the portion of the wavefunction 
containing the terms A e , which we can write in the form: 

&l(«i,*a) = e 2klZl+tk2Z2 + ae~ iklZl+ik2Z2 (59) 

where the energy is equal to E = k\ + k 2 and could 
be complex. Similar to the single-particle solutions, the 



presence of the imaginary part in the energy reflects the 
fact that the two-particle state stays a finite amount of 
time inside the system. Applying boundary conditions at 
z = and z = d subsequently generates four equations 
relating a, /3,7 where one of them is redundant. Their 
solution reduces the wavefunction to 

<f>A{zi,Z2) = e iklZl+ik2Z2 J- h±Ao-ikiZi+ik 2 z 2 



ki - 1 



(60) 

+ik\z\ — ik 2 z 2 



kx-1 

where 7 = k2 ~^\ e 2lk2d ■ A similar expression results for 
the portion of (f>(zi, z 2 ) containing the B e terms, once the 
boundary conditions at z — and z = d are applied: 



4>b{zi,z 2 ) = - l - 



2 Zi+ik\Z2 _|_ ^jL_t_!-p~ ik 2 zi+ikiz 2 

k2 — 1 



I k 2 4- 1 ^— ifc 2 ^ 1 — ifc lZ2 ^+ik 2 zi—ikiz 2 



where 7' = j^je 2lkld , and t is a coefficient to be de- 
termined from the boundary condition at Z\ = z 2 . To 
find t, it is convenient to re-write each of ther terms in 
$>A,B as a product of relative coordinate (r = z 2 — z\) and 
center-of-mass coordinate (R = (z\ + z 2 )/2) functions, 



4> A (R,r) 



^ipR—iqr _j- ^ c ~ iqR+ipr 



ki - 1 



(61) 



u 1 1 

j 1 Z c -ipR+iqr 1 iqR-ipr^ /g2) 

ki — 1 



MR,r) = ( e ^+^ + 

t «2 — 1 



(63) 



7 



1 



e -ipR-iqr _|_ y & -iqR-ipr\ 



where p — (k% + k 2 ) and q = (k\ — k 2 )/2. The bound- 
ary condition at z\ — z 2 leaves the center-of-mass parts 
of the wavefunction unaffected, but yields the following 
condition on the relative coordinates, 

d r (t>(R,r)\ r=0 + = —(f>(R,r)\ r =o+- (65) 

where <fr — 4>a + 08 is the total wavefunction in the tri- 
angular region. We should satisfy this boundary condi- 
tion separately for each of the center-of-mass momentum 
terms e pR , e ±zqR in the total wavefunction. This leads 
to three independent equations (one out of four is redun- 
dant). However, we introduce a new parameter (t 1 ) to 
simplify the equations, which turns them into four equa- 
tions: 



t = 



ki 



IK 



IK 



ki 
it' 

t' 



k 2 
ki- 



ki + k 2 



ki-1 



J2ik\d 



_ £ e 2ik 2 d 



(66) 
(67) 
(68) 
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which can be written in the following short form: 

c 2ikjd = ( k i + !) 2 -Q ( fc » ~ k j + iK )( k i + k 3 + iK ) (gg) 

{J\% 1) . , (A:? fcj itz^i^ki ~\~ kj zk) 

where «, j can be (1,2). These are transcendental equa- 
tions for (ki,k 2 ), which generate the spectrum of two 
interacting particles. We can also write the wave func- 
tions (0 = 4>a(zi,z 2 ) + 4>b(zi,z 2 )) in the region (0 < 
z\ < z 2 < d) in a more compact way, by using the single 
particle solutions T]k(z) = sin(fcz) + ikcos(kz): 

4 e ik 2 d 

<\>a{z\,z 2 ) = — —riki{zi)rik 2 {d - z 2 ) (70) 

k\ — 1 fca + 1 

</>b(zi,z 2 ) = -, --, — — rj k2 (zi)7j kl (d ~ z 2 ). (71) 

k 2 — 1 k\ + 1 

It is interesting to note that in the limit of strong in- 
teraction (either for positive or negative k), the solutions 
are very similar to the non-interacting case. The reason 
can be seen from the transcendental Eqs. in the 

limit n — ► ±oo. We then recover the same characteris- 
tic equations e 2lkd = (j^j) 2 for both wavevectors as the 
non-interacting case, Eq. ([oTj). We should note that there 
are some trivial solutions to the transcendental Eqs. (|69|) . 
which do not have any physical significance. For exam- 
ple, equal wave vectors k\ — k 2 . Although one can find 
such wave vectors, this solution is readily not a solution 
to Eq. ([53|) . since it does not satisfy the interacting part 
(this solution only contains center of mass motion). One 
can also plug back the wave vectors into wave function 
and arrive at a wave function equal to zero everywhere. 
Another example is when one of the wave vectors is zero. 
In this case, one can also show that the wave function 
is zero everywhere. If (ki,k 2 ) are solutions to the tran- 
scendental equations, then (±ki,±k 2 ) are also solutions 
with equal energies. In next two sections, we investigate 
non-trivial solutions to the transcendental equation for 
two particles and discuss the related physics. 

C. Solutions close to non-interacting case 

The transcendental equations allow a set of solu- 
tions with the wavevectors close to two different non- 
interacting modes say (m, n). In the non- interacting 
regime, any mode can be populated by an arbitrary num- 
ber of photons. However, once the interaction is present, 
photons can not occupy the same mode and therefore, 
the photons will reorganize themselves and each acquire 
different modes. Fig. [5] shows a normal mode wavefunc- 
tion of a non-driven system in both non-interacting and 
strongly interacting regime (nd 1). The wave func- 
tion has a cusp on its diagonal and diagonal elements are 
depleted for both repulsive and attractive strong interac- 
tion. 



This is a manifestation of fermionization of bosons in 
one dimensional system in the presence of strong interac- 
tion [H, [33] • Such solutions can exist both for repulsive 
and attractive interactions. However, we note in the case 
of attractive interaction such solutions are not the ground 
state of the system and solutions with lower energies ex- 
ist which will be discussed below. We later argue that 
indeed on the repulsive side, the anti-bunching behavior 
of a driven system is due to the repulsion of the photons 
inside the medium. We can also estimate the energy of 
such modes which is always positive. In the strong inter- 
acting regime, particles avoid each other and therefore, 
their energy of a strongly two interacting bosons E(m, n) 
will be equal to the energy of a system which has two 
non-interacting bosons, one in state m and the other in 
state n. This is shown in Fig. [5J where by increasing 
the interaction strength the energy of interacting parti- 
cles reaches that of the non-interacting particles. As we 
pointed out in the previous section (Scc lVI A)) , the energy 
of modes {E{n)) in an open box has an imaginary part 
which represents how fast the particle leave the system. 
However, for large systems (d ^> 1), this decay is very 
small compared to the energy of the mode and one can 
approximate the energy of an open system by that of a 
closed box (i.e. E(n) ~ (^f) 2 )- Therefore, the energy of 
two strongly interacting photons (nd 3> 1), in the limit 
of large system (d^> 1), will be given by: 

E(n,m)c(-) +(-) (72) 

We note that our strongly interacting system is char- 
acterized by the parameter nd which is the same 7- 
paramctcr conventionally used for interacting ID Bose 
gas. More precisely, the 7-parameter which is the ratio 
of the interaction to kinetic energy can be simplified in 
our case for two particles: mud/2 — nd/4. 

-g Kd-50 

Lr\ r% L'\ 

FIG. 8: The amplitude of the two-photon wavefunction for 
(m=l, n=2) mode when the system is not driven. By increas- 
ing the interaction, photons self-organize inside the medium 
and exhibit anti-bunching (depletion of diagonal elements). 
For this plot: d = 30. 



D. Bound States Solution 

For attractive interaction (k < 0), the mode equa- 
tion ([55)1 admits solutions which take the form of pho- 
tonic bound states. Specifically, in the reference frame of 
the center of mass, two particles experience an attractive 
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FIG. 9: Energy of two-photon states: by increasing the in- 
teraction strength (nd 1) the energy of interacting par- 
ticles (red:E(l,2) and blue:E(l,3)) reaches that of the non- 
interacting particles (black). For large system (in this case 
d = 30 ^> 1), the energy limit is equal to energy of particles 
in a closed box (green). The error bars show that imaginary 
part of the energies. 



delta function interaction —2\k\S(z 2 — zi) — > — v2|«|£(r), 
which allows one bound state in the relative coordinate. 
Therefore, the part of the wavefunction describing the 
relative coordinate roughly takes the form e I<? ' r ', where 
the relative momentum q = (ki — k 2 )/V2 — i\n\/\/2 is 
imaginary and its energy is about — k 2 /2. On the other 
hand, the center of mass momentum can take a discrete 
set of values that are determined by the system boundary 
conditions. We find that the center of mass solutions can 
be approximately described by two different types. The 
first type is where the real part of each photon wavevector 
roughly takes values allowed for a single particle in a box, 
such that ki~ (^)+i§ and fc 2 ^ ( I r)-*f- In this case 
the center of mass has wavevector p = fci + k 2 — 2(^| L ). 
The corresponding energies for these states are 



Et 



/mr\ 2 n 2 

2( t) -y 



(73) 



Here, the first term on the right corresponds to the en- 
ergy of the center of mass motion, and the second term 
corresponds to the bound-state energy of the relative mo- 
tion. Fig. [TO] shows that the energies estimated in this 
way agree very well with the exact values obtained by 
solving the transcendental Eqs. ([uT)]) . 

The second type of solution allowed for the center of 
mass motion is where its energy approximately takes a 
single-particle value, 2 hm) = ) 2 wnere p = &i + &2 — 
y/2(HS). Therefore, the momentum of individual parti- 
cles will be given by k\ ~ (7^) + if > fc 2 ^ (7^) ~ *§ 
and the energy of this paired composite can be estimated 
as 



Et 



nir\ 2 k 2 



d ) 



(74) 



Again, the estimated energies agree well with exact so- 
lutions, as shown in Fig. [T2J We note that some of the 



FIG. 10: Energy of bound states versus strength of nonlin- 
earity. Green (solid) curves are obtained by solving transcen- 
dental Eqs. (|69[) . Blue (dashed) curves are estimated based 
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FIG. 11: The amplitude of the two-photon wavefunction for 
a bound state when the system is not driven. By increas- 
ing interaction, photons becomes more bunched. (^1,^2) — 
(IHL) ± j«. For this plot: d = 60, n = 3. 



estimated allowed energies for the two types of center 
of mass solutions coincide (e.g., the lowest lying energy 
level in Fig. [TO] and Fig. [IT 



The energies of this series of bound states decrease with 
increasing strength of nonlinearity Now, suppose we 
drive the system with a coherent field of fixed frequency 
(5, while varying n. The system is expected to display a 
set of resonances as |«| is increased, each time S is equal 
to some particular bound state energy E h n . This effect in 
fact gives rise to oscillatory behavior in the correlation 
functions as a function of k, as we will see later (Fig. 
QUa) and (b)). 

The wavefunction amplitude of a typical bound state 
is shown in Fig. 1131 Due to the attractive interaction, 
diagonal elements z\ = z 2 become more prominent as |k| 
increases, indicating a stronger bunching effect for the 
photons, and these states become more tightly bound in 
the relative coordinate. The center of mass of the bound 
states can acquire a free momentum that is quantized 
due to the system boundary conditions (e.g., k ~ nir/d). 
Fig. [TO] shows the wavefunction of the third bound state 
(n=3). The three peaks evident for large |«| reflect the 
quantum number of the center of mass motion. 
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two particles. In the limit of hardcore bosons where k — > 
oo, the expressions can be simplified since t, t! = — 1 and 
e lkd = j^j for both k = k\^. Then, the two components 
of the wavefunction 4>a and 0e take very similar forms, 



(fcl - l)(fc 2 - 1) 

-4 

(k x - l)(k 2 - 1) 



Vk 1 {zi)vk 2 (z2) (79) 
Vk 2 (zi)Vk 1 (z 2 ). (80) 



The generalization to the many-body solution is straight- 
forward for the hardcore boson case (also see Ref. [381]): 



FIG. 12: Energy of bound states versus strength of nonlin- 
earity. Green (solid) curves are obtained by solving transcen- 
dental Eqs. (|69[) . Blue (dashed) curve are estimated based 
on E b n ~ (if) 2 - ^. In this plot d = 30. 
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FIG. 13: The amplitude of the two-photon wavefunction for 
a bound state for a non-driven system. By increasing interac- 
tion, photons becomes more bunched, (ki, fo) — ■ 
For this plot: d — 60, n — 3. 



E. Many-body problem 

In this section, we obtain the general solution for 
the many-body case. For the many-body system, the 
Schrodinger equation takes the form 



1 d 2 

E<f>(zi,...,z N ) = -— 22-q -s<f>(zi,.-,z N ) (75) 

i 1 

+ J~] 2K(f>(zi > ...,z N )S(zj - Zj),(76) 

<i,j> 



[Zl,Z 2 , 



,z N ) 



det Vj(zk) 

l<j,k<N 



(81) 

Similar to two-body solution, we note that such solutions 
are present both for positive and negative k. Since the 
system is one dimensional, strong interaction leads to 
fermionization of bosons ( in this case photons) [H, [3?J • 
We can also extend the many solution for an arbitrary 
interaction strength, following Refs. (l8l. |39| . Similar to 
the two-body case, we can construct the general many- 
body wave function of the form: 



(zi,z 2 , . 



.z N 



E 

p 



(82) 



where the first sum is over forward and backward going 
waves (e = ±1) and the second sum is over different mo- 
mentum permutations of the set {k} — (k±, k 2l fc^r), 
therefore there are 2^^! terms. We can find Bp co- 
efficient by requiring ^ p Bp Yl%<j e lep i kp i Zi to be solu- 
tion to the Schrodinger equation (Eq 175)1 . We can write 
these coefficients in a compact way according to Gaudin 

m & Bp = (l + 6pA< i\ fc J with the total 

energy E = Yli^i- Now, we apply the boundary con- 
dition Eq. ([77jl which relates coefficient A e . For a given 
momentum permutation P = (pi,P2,, ---tPn), by consid- 
ering the terms corresponding to different signs of e Pi , 
the boundary condition requires A e to satisfy equations 
of the form 



where < i,j > indicates pairs of particles. The open 
boundary conditions for the many-body problem are 
given by 



(zi, z N ) - i—<f>(zi, z N ) 

OZi 

d 

>(zi,...,z N ) +i — (j)(z 1 ,..,z N ) 



Zi=0 



(77) 



0. (78) 



Before presenting the general many-body solution, we 
first study the limit of very large interaction strength for 



(1 + £p; kp t )A ei t ,, e> 

(1 - fpAi M^, ..(-£„,). 



n 

« n 



IK 



£pi kpi 



k Pi 



= 



The above equations can be satisfied by the following 
solution for A f 



A, 



n 

i<3 



IK 



tiki 



N 

<ikj n 

J ■> ' 771 — 1 



6 m k r . 



(83) 
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Therefore, the wavefunction can be written as: 

<j> (zi,Z2, ...,zn) = (84) 
N / 1 \ 

J^J I 1 — ] e l l ( e Pl k Pl x l+--+ € PN k PN XN ^ 



n 



P m=l 
1 



tiki ~\- €jkj 



Similar to two-body case, we have to subject this solution 
to the boundary condition at other end (i.e., z = d) to 
determine the momenta fcj's. This condition yields the 
transcendental equations for momenta: 



k,d _ (ki + I) 2 TT (^» _ + + % + i/c) 



— TT 

n 2 11 



(fcj l) 2 . .^ (fcj kj in)(k% + fej 



(85) 



If we assume only two particles in the system, one can 
easily verify the the above transcendental equations re- 
duce to two-body transcendental equation derived in the 
previous section (Eg. ([55])). 



VII. QUANTUM TRANSPORT PROPERTIES 

In this section, we investigate transport properties of 
the photonic nonlinear one-dimensional system in the 
regimes of attractive, repulsive, and absorptive interac- 
tions between photons. We present numerical solutions 
for the transport of photons incident from one end of 
the waveguide (a driven system), while using the analyt- 
ical solutions of the non-driven system (Sec. IVIj) to elu- 
cidate the various behaviors that emerge in the different 
regimes. 



A. Repulsive Interaction (« > 0) 



field. In Fig. [T4] the transmission of the single-photon 
intensity 



Tx = 



(l|*t(0)* + (Q)|l)' 



the transmission of the two-photon intensity 
(2|*t(d)*+(d)|2) 



To 



(2|*t(0)*+(0)|2) 



(86) 



(87) 



and the transmitted correlation function gi(z = d, r = 
0) is shown as the system evolves in time. The system 
reaches its steady state after a time of the order of the 
inverse of the system bandwidth (Sec. IIII[) . In fact, T\ 
coincides with the linear transmission coefficient of the 
system in the absence of the nonlinearity. 

(r iD ,A r A 2 )/r=(0.1, 1,4.1), OD=160 
0.7 ■ n O.2 




2 T(;t/d) 3 3 



0.15 

o.i ■; 

0.05 Y 
,0 



FIG. 14: g-iir = 0) reaches the steady-state after a time inter- 
val which is set by the bandwidth of the system, one-photon 
state (green) is partially transmitted while the transmission 
of the two-photon state (red) is further suppressed due the 
nonlinear dispersion. This has been generated for a system 
with Tid/Y = 10% an OD=160. 



We first study the quantum transport properties of the 
system in the dispersive regime where the nonlinearity 
coefficient is almost real and positive (k > 0), such that 
photons effectively repel each other inside the system. 

We assume that a weak coherent field is incident to 
the waveguide at one end, z = 0, with no input at the 
other end, z — d [similar to Fig. HJb)]. We fix the de- 
tuning of the input field to <5o = (ir/d) 2 , which corre- 
sponds to the first transmission resonance in the linear 
regime (Sec. IIII[) . Because we have assumed a weak in- 
put field, we can apply the techniques described in Sec.lVl 
to describe the transport. Our numerical techniques for 
solving these equations are given in Appendix [Bj While 
the numerical results presented in this and the follow- 
ing sections are evaluated for a specific set of parameters 
(system size, detuning, etc.), the conclusions are quite 
general. Numerically, we begin with no photons inside 
the medium, and evaluate quantities such as the trans- 
mission intensity and correlation functions only after the 
system reaches steady state in presence of the driving 



First, we note that the single-photon wave function is 
not affected by the presence of the nonlinearity and will 
be perfectly transmitted in the absence of linear losses. 
Thus, in our truncated Hilbert space, the only subspace 
affected by k is the two-photon wave function, which is 
shown in Fig. [15] We clearly observe that the nonlin- 
earity causes repulsion between two photons inside the 
system, as the wave function along the diagonal Z\ = Z2 
becomes suppressed while the off-diagonal amplitudes be- 
come peaked (indicating the de-localization of the pho- 
tons). This behavior closely resembles that of the natural 
modes of the system, as calculated in Sec. I VII A simi- 
lar behavior involving the "self-organization" of photons 
in an NLSE system in equilibrium has been discussed in 
Ref. [|. 

In the presence of linear absorption (discussed in 
Sec. IIII[) . the system will not be perfectly transmitting 
even on resonance, and therefore in a realistic situation 
the transmittivity will be less than one (T\ < 1). Note, 
however, that such absorption would result in a classical 
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FIG. 15: Two-photon wave function |0(zi,Z2)| exhibiting de- 
localization. We have assumed no dissipation (r' = V — 0) in 
this plot, d — 30 for different values of k. 



output given a classical input. Significantly, in the pres- 
ence of a nonlinearity, we find that the output light can 
acquire non-classical character. Specifically, the trans- 
mitted light exhibits anti-bunching (g2(z = d, t = 0) < 
1), which becomes more pronounced with increasing nd 2 
(Fig. HBjl . This effect partly arises from the suppression 
of transmission of two-photon components, due to an ex- 
tra nonlinear phase shift that shifts these components 
out of transmission resonance. In fact, these components 
are more likely to get reflected, which causes the reflected 
field to subsequently exhibit bunching behavior. We note 
that this effect is similar to photon blockade in a cavity 
(e.g., see Refe.ja El, In addition, additional anti- 

bunching occurs due to the fact that two-photon compo- 
nents inside the system tend to get repelled from each 
other. This effect arises due to the spatial degrees of 
freedom present in the system, which is fundamentally 
different than switching schemes proposed in optical cav- 
ities (e.g., Refs.Jl, Jl, 03) or wave-guides coupled to a 
point-like emitter [3|, Il7| . In the limit where k — > oo, 
the transmitted field approaches perfect anti-bunching, 
ff2 (rf,r = 0) = 0. 
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FIG. 16: gi(r — 0) as a function of nonlinearity. For large 
system sizes d ^> 1, the anti-bunching of the system scales 
with Kd 2 . 

In an experimental realization, the requirement to see 
the photon repulsion (nd 2 > 40) for a system with 



Tid/T = 10%, would be a coherent optical length of 
d ~ 40 when A2/I 1 = 1. Therefore, at least an optical 
density of OD ~ 160 is needed for T x ~ 20%. The anti- 
bunching in the transmitted light is more pronounced as 
the optical density increases, which increases the effective 
system finesse (Fig. ITT1) . 




OD 

FIG. 17: Repulsive photons: Correlation function gi (r = 0) 
of the transmitted light when the frequency is set to the single- 
photon transmission resonance with Ti ~ 20% and = 5. 



B. Attractive Interaction (k < 0) 

In this section, we study the quantum transport prop- 
erties of the system in the presence of dispersive non- 
linearity with negative coefficient. Contrary to the semi- 
classical prediction, we show that the second-order corre- 
lation function of the transmitted field oscillates as func- 
tion of nonlinear interaction strength and can exhibit 
both bunching and anti-bunching. We explain the ori- 
gin of this behavior in terms of the analytical solutions 
obtained in Sec. IVIB1 

In Fig. UST a). we plot 32 (t = 0) for the transmit- 
ted field versus nd. Initially, the system exhibits anti- 
bunching behavior for small values of \n\d which indi- 
cates that multi-photon components tend to switch them- 
selves out of transmission resonance. However, as we 
increase \n\d, oscillations develop in the correlation func- 
tion, exhibiting strong bunching behavior at particular 
values of Kd. Thus, unlike the repulsive case, a com- 
peting behavior arises between the photon switching ef- 
fect and the resonant excitation of specific bound states 
within the system, as we describe below. In particular, 
the bound state energies ~E\ decrease quadratically with 
changing K, according to Eq. ([74")) or Eq. ([75]) . which is 
shown in Fig. [TST b). For a fixed detuning 6, the oscilla- 
tion peaks (where gi is largest) correspond to situations 
where the energy of a bound state becomes equal to the 
energy of two incoming photons [E\ = 28). This effect 
is further confirmed by examining the two-photon wave 
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function at each of these oscillation peaks (Fig. [LSk). We 
clearly observe that these wave functions correspond to 
the bound states calculated in Sec. IVIBl Similar to Fig. 
[HI and Fig. [T3J it is readily seen that the wave functions 
at these peaks are localized along the diagonal, indicating 
a bound state in the relative coordinates and leading to 
the bunching effect in transmission. On the other hand, 
an increasing number of nodes and anti-nodes develop 
along the diagonal for increasing which are associ- 
ated with the higher momenta of the center-of-mass mo- 
tion. We note that such resonances deviate significantly 
from the the semiclassical picture, where anti-bunching 
was predicted for both positive and negative nonlinear- 
ity. We also note that in cavity QED systems this effect 
is not present since these systems are single-mode. 




FIG. 18: (a) Output correlation function g^ir — 0) as a func- 
tion of nonlinearity: When the negative nonlinear strength is 
changed to higher values and t/2 exhibits resonances at certain 
values of Kd ~ (0, 6, 10, 14, ...). In this plot the system size is 
d = 30, however, for other system size same behaviors were 
observed around similar values of kcL. The two-photon wave- 
function (|0(2i,2;2)|) for four values of nonlinearity is shown, 
(b) Corresponding bound state energies (green-solid) which 
become resonant with incoming photon energy (black-dotted) 
for specific nonlinearities. We have assumed no dissipation 
(r' = r = 0) in these plots. 



The experimental requirement to see such behaviors 
is more stringent than the photon repulsion in the pre- 
vious section. For example, if we want to observe the 
second photonic bound state (nd > 5) for a system with 
Tid/T = 0.2, the coherent optical length should be at 
least d ~ 200 when /S.2/T = — 5. To achieve a reasonable 
signal (linear transmission T\ = 1%) an optical density 
of OD = 3500 is needed. Importantly, however, we have 
shown that the presence of bound states inside the non- 



linear medium can be probed with classical light, simply 
by examining higher-order correlation functions in the 
output field, rather than sending in complicated quan- 
tum inputs. 



C. Dissipative Regime (n — i\tt\) 

In this section, we study the transport properties of 
the system in the presence of nonlinear absorption, and 
calculate its effect on the transmitted light and its corre- 
lation functions. 

A purely absorptive nonlinearity arises when the de- 
tuning A2 is set to zero in our atomic system (see 
Fig.nja)). This nonlinear loss also leads to anti-bunching 
in the transmitted field, as multi-photon components be- 
come less likely to pass through the waveguide without 
being absorbed. Linear absorption, on the other hand, 
affects transmission of single- and multi-photon compo- 
nents equally. Fig. [19] and Fig. [20l show how two-photon 
and one-photon states are transported differently in the 
nonlinear absorptive system (realistic linear losses are in- 
cluded in this calculation). 



(r, r ,,A,,A, 1 )/ T=(0.1,0 ) 2-7), OD=70 
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FIG. 19: giir = 0) reaches the steady-state after a time 
interval which is set by the bandwidth of the system, one- 
photon state (green) is partially transmitted while the two- 
photon state (red) is strongly attenuated due the nonlinear 
absorption. This plot has been generated for a system with 
riu/r = 10% an OD=70. 



We note that the two-photon wavefunction is atten- 
uated due the nonlinear absorption, while it is not de- 
formed, as shown in Fig. 1211 In an experimental real- 
ization of such a system with Tid/T = 10%, an optical 
coherent length of d ~ 20 is enough to yield a relatively 
large anti-bunching (g 2 < 0.3). In order to have high 
transmission (T\ = 20%) for single photons an optical 
density of OD ~ 70 is required. 

All of the physics related to the photon correlation 
function is described again by product of the coherent op- 
tical length and the nonlinearity coefficient (since 
the nonlinear absorption is equal to the nonlinear absorp- 
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This technique allows us to study the system even in 
regimes where nonlinearities are significant even at a few- 




200 



200 



FIG. 20: (a) one-photon state is partially transmitted (Ti) 
while the two-photon state transmission (T2) is suppressed 
due the nonlinear absorption. This suppression is more pro- 
nounced for higher optical density and cooperativity. (b) Cor- 
relation function gi(r = 0) of the transmitted light when the 
frequency is set to the single-photon transmission resonance 
with T ~ 20% and ^ = 0. 



tion coefficient times the length of the medium). How- 
ever, for a fixed optical density, since the nonlinear tran- 
sition is on resonance, the magnitude of the nonlinear 
coefficient \k\ is enhanced compared to the nonlinear dis- 
persive case. We note that in the presence of nonlinear 
absorption, one has to also consider the effect of accom- 
panied noise. However, the effect of noise for an ensemble 
of many atoms which are driven by a weak laser field, is 
negligible, and therefore, using the NLSE with a decay 
term is sufficient and consistent. A rigorous demonstra- 
tion of the validity of such approximation is the subject 
of further research. 



VIII. CONCLUSIONS 

We have developed a technique to study few-photon 
quantum dynamics inside ID nonlinear photonic system. 




FIG. 21: In the presence of the nonlinear absorption = 
0), the two-photon wave function is strongly suppressed com- 
paring to the absence of the nonlinear absorption (A2 3> T). 
These plots has been generated for a system with Tid/T = 
10%, Tj ~ 20% and OD=70. 



photon level, where we find that the behavior of the sys- 
tem deviates significantly from estimates based on clas- 
sical formalism. Specifically, when the system is driven 
by classical light, the strong optical nonlincarity man- 
ifests itself in the correlation functions of the outgoing 
transmitted light. In particular, when the interaction 
between photons is effectively repulsive, the suppression 
of multi-photon components results in anti-bunching of 
the transmitted field and the system acts as a single- 
photon switch. In the case of attractive interaction, the 
system can exhibit either anti-bunching or bunching, as- 
sociated with the resonant excitation of bound states of 
photons by the input field. These effects can be observed 
by probing statistics of photons transmitted through the 
nonlinear fiber. 
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APPENDIX A: EIT-BANDGAP 

In this appendix, we show that how in an EIT system, 
where the control field is a standing wave, a band gap 
structure can be developed. In particular, we show the 
presence of transmission resonances at the band gap edge 
by taking into account the full expression for the atomic 
susceptibilities. We show that at the band gap edge, we 
recover that same resonances that we presented in the 
main text for small detunings. 

We consider a A-level scheme, where a standing control 
field has coupled the forward- and backward-going probe 
together, similar to Fig[T] without the nonlinear transi- 
tion (c-d). Following [421] . we assume the noises to be 
negligible, and therefore, the atomic equations of motion 
to the leading order in £± are 

d t°t b = - T/2)a+ b + iCla ac + igV2^£+(Al) 

d tKb = +(iAi-r/2)a- 6 + iOa oc + ipV2^£-(A2) 

d t (Tac = -70<7 ac + i^ab + ^Kb ( A3 ) 

and the evolution equation of the photonic fields are writ- 
ten as: 

(d t + cd z )£+ = iAK£+ + igV2^n a+ b (A4) 
(d t - cd z )£- = iAK£- + igV^noa^. (A5) 

The wave- vector mismatch can be ignored by including 
a small shift in the two-photon detuning. By taking the 
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Fourier transform of the atomic equation of motion, one 
can solve for atomic polarization and obtain the self- and 
cross-susceptibilities. We can define a unit length based 
on the absorption length L a t, s = 97r g^„ o , and write the 
field equation as: 



A/f=-50, OD=5, n=20 



d~ z £+ = iA 3 £+ + i Xs (6)£+ + i X c(S)£- (A6) 
-d 2 £„ = iA 3 £_+i Xs (S)£-+i X c{S)£ + 

where the self- and cross susceptibilities and the detuning 
are given by: 



Xs(8) 
Xo(S) 
A 3 



r rT 







r> 2 



T' TTo + 2fL 2 

.r r> 2 



2n 2 

s r 

f r 2ng 2 n ~ 



V T'Tq 

A, r 2 



(A7) 
(A8) 
(A9) 



where V = T/2- iAi - iA 3 , T = 70 - iA 3 and A 3 
is the two-photon detuning of the probe from the pump 
field which is related to the dimensionless two-photon 
detuning in the main text (A3 = 2^5 < Aj,). We 
note that in most cases, A3 is very small for slow group 
velocities (tt^s — = (?t) 2 ^t — <C 1), and therefore the 
corresponding term can be neglected for simplicity. 

In order to obtain transmission and reflection co- 
efficient, one should solve the couple mode equations 
Eq. (|A6p with proper boundary conditions. Therefore, we 
consider a system which is driven with a weak coherent 
field at [z = 0). Therefore, the boundary conditions can 
be set to, 



£+(z = 0) = £ a 
£_{z = d) = 0. 



(A10) 
(AH) 



We evaluate the transmission coefficient (£+(z — 
d)/£o), and the reflection coefficient (£-(z = 0)/£q) by 
numerical methods using BVP5C in Matlab. In partic- 
ular, we are interested in the Raman regime, in other 
words the detuning is very large |Ai| ^> T and also we 
assume Ai < 0. First, we consider the case where the 
EIT width is smaller than the one-photon detuning, i.e. 
fl -C |Ai|. Fig |2"21 shows the reflectivity and transmit- 
tivity of the system for different optical densities. In 

the regime with low optical density, the spectrum corre- 
ct 2 

sponds to a shifted Raman transition at A3 ~ 2^- and 
an EIT window around A3 ~ — Ai . In higher optical den- 
sities, the system develops a band gap for — < A3 < 0. 
Figf22] shows that in media with higher optical density, 
the band gap becomes more prominent. 

As we discussed in the main text, we are interested 
in the band gap edge where the transmission peaks 
are present and the system acts like an effective cavity. 
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FIG. 22: By increasing OD, the band gap structure becomes 
more pronounced 



Figl23l shows a close-up of the transmittivity and reflec- 
tivity spectrum in Fig|2"2Tc) at the band gap edge. We 
can observe that several resonances occur at the edge due 
to the finite size of the system. By positioning at the one 
of the transmission peaks, the system behaves as an ef- 
fective cavity, where the decay rate of the cavity will be 
given by the width of the transmission peak. Therefore, 
the present results, including the full susceptibilities of 
the system, is consistent with the model presented ear- 
lier where we had approximated the system to be around 
A 3 = 0. 

We add that alternatively, one can assume a strong 
control field so that the EIT windows would be smaller 
than the one-photon detuning fl > |Ai|. Similar to the 
previous case fi < |Ai|, the system develops a band gap. 
As shown in Fig|24l the band gap is formed between 
— £1 < A3 < fi, similar to modulated EIT with AC-stark 
shift as discussed in Ref . [241] . 



APPENDIX B: NUMERICAL METHODS 

In this section, we describe the numerical methods that 
have been used to simulate the evolution of the photonic 
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FIG. 23: Transmission resonances at the edge of the band gap 



A,/T=-20, OD=1000, £2=100 




200 



FIG. 24: Band gap structure for strong control fields (fi > 
Ai) 



quantum state and the related correlation functions, in 
the limit where we truncate the Hilbert space to two 
photons or less. The partial differential equations for the 
one-photon and two-photon wave functions (|27l [28)) are 
turned into difference equations by discretizing space and 
time, and are evolved forward in time using the Du Fort- 
Frankel scheme [43| . This algorithm is is explicit in time - 
i.e., the next time step function is explicitly given by the 
past time function - and is also unconditionally stable. 
We note that the system under investigation is open and 
it is driven out of equilibrium, therefore, conventional 
analytical methods for approaching the NLSE such as 
Bethe ansatz or quantum inverse scattering [1 91 ] are not 
applicable here. 

The one-photon wave function can be easily integrated 
and solved analytically. However, we describe how to ob- 
tain the one-photon wave function numerically and then 
generalize this technique to obtain the two-photon wave 
function. First, we mesh space and time and reduce 
the differential equations to a difference equation. If we 
choose the time step k and the space step h, the dis- 
cretized time and space will be x — z/h and s = t/k 
and the system length d = Nh. Then following the Du 



Fort-Frankel scheme [43(, the evolution equation takes 
the form: 



(x, s + 1) - 0(x,s - 1) 
2k 



— Ms + M) +0(3-1, a ) 
d(x,s + l)-9(x,s-l)] (Bl) 



where the position take all values inside the boundary 
(2 < x < N — 1). By rearranging the above equation, the 
explicit form of the equation can be obtained 



( 1 + J ^)0(x,s + l)=9(x,s-l) 

' k [6[x + l,t) + 0(x-l,t)-d(x,s-l 



(B2) 



mh 2 



Therefore inside the boundaries, the wave function at 
time s + 1 can be obtained knowing the wave function at 
time s and s — 1. The boundary condition at z — -i.e. 
x — 1, will be given by 



9(1, 
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6(2, 
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.9(2, t- 
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(B3) 



2mh 



Or equivalcntly, 



61(1, 8+1) 



"+(-1+2^)^(2^ + 1) 
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and similarly for the boundary condition at z = d -i.e. 
x = N, we have 



9(N,s + l) 



H + sbW^-M + i) 



2 mh 
1 
2 



(B5) 
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Therefore, by having the above boundary conditions and 
the initial condition 9(x,s = 1) = 0, the wave function 
can be calculated at any time inside the boundaries (2 < 
x < N—l). The order of accuracy of the Du Fort-Frankel 
scheme is given by 0(h 2 ) + Q(k 2 ) + 0(k 2 h~ 2 ) and it is 
consistent as k/h tends to zero [43]. 

Similarly, we can write a difference equation for the 
two-photon wave function. The (5-interaction can be ap- 
proximated by a Gaussian distribution. The space do- 
main is meshed so that A.z% = Az 2 = h. The evolution 
equation for the two-photon wave function reads 



1 



2ik 
mh 2 



(x,y, s + l) 



(B6) 
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The boundary condition at z = will be given by 

f%, S +l) = *(1,*M + 1) + *(2,IM + 1) (B7) 

_ , 0(2, + 
2to/i 

where cr is the length scale characterizing the distance 
of the two-photon interaction. Approximating the delta- 
function with a Gaussian is valid if tr rf. On the other 
hand, we should have h <C o so that the Gaussian func- 
tion would be smooth. Or equivalently, 

, (1>p> , + 1) = M^ + l) + H^(2, + 

2 2rad 

(B8) 

and similarly for the boundary condition at z = d, we 
have 

Q = </>(N,y,s + l)+<l>(N-l,y, S + l) 



+ i 
which gives 

4>(N,y,s+l) 



</>(N,y,s + l)-<l>(N-l,y, S + l) 
2mh 



- + — 

2 ~ 2mh 



(BIO) 



Once the wave function is known at any point in time and 
space, we can evaluate the correlation functions. In par- 
ticular, the two-photon correlation function gi{d,T = 0) 
is given by Eq. (|40p . where the first and the second deriva- 
tives at anytime are given by the following expressions, 



d (1) (f>{d, d) 
5 (1) 5 (2) 0(d, d) 



2mh 
1 



[<fi(N,N)-ct>(N-l,N)](Bll) 
:[<j>(N,N)-<f>(N-l,N) 



Am 2 h 2 

- <f>{N,N-l) + c/>(N-l,N-l)]. 



Note that in evaluation of g% (r), once the first photon 
is detected the two-photon wave function collapses to 
zero. This seems to be contradictory with the driven 
boundary condition Eq. (|55|) where the two-photon state 
at the boundaries is proportional to the one-photon wave 
function which is not zero. This apparent inconsistency 
occurs because we have neglected higher number pho- 
ton states in our truncation. However, this inconsistency 
only leads to higher order corrections to gi (r) in the input 
field amplitude a, which is assumed to be weak (a<Cl). 



